Load libraries need to process the data
# load library
library(raster)
## Loading required package: sp
library(rhdf5)
library(rgdal)
## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-3
# check working directory
getwd()
## [1] "/Users/cassondrawalker/Documents/data/NEONDI_2016/monday_pm"
Open data file
# Open file object by defining file path
# "tab" key allows you to select an element in a folder
# ../up one folder
f <- "../NEONdata/D17-California/TEAK/2013/spectrometer/reflectance/Subset3NIS1_20130614_100459_atmcor.h5"
# View h5 structure
h5ls(f)
## group name otype dclass
## 0 / ATCOR_Input_File H5I_DATASET STRING
## 1 / ATCOR_Processing_Log H5I_DATASET STRING
## 2 / Aerosol Optical Depth H5I_DATASET INTEGER
## 3 / Aspect H5I_DATASET FLOAT
## 4 / Cast Shadow H5I_DATASET INTEGER
## 5 / Dark Dense Vegetation Classification H5I_DATASET INTEGER
## 6 / Haze-Cloud-Water Map H5I_DATASET INTEGER
## 7 / Illumination Factor H5I_DATASET INTEGER
## 8 / Path Length H5I_DATASET FLOAT
## 9 / Processing Version H5I_DATASET STRING
## 10 / Reflectance H5I_DATASET INTEGER
## 11 / Shadow_Processing_Log H5I_DATASET STRING
## 12 / Sky View Factor H5I_DATASET INTEGER
## 13 / Skyview_Processing_Log H5I_DATASET STRING
## 14 / Slope H5I_DATASET FLOAT
## 15 / Slope_Aspect_Processing_Log H5I_DATASET STRING
## 16 / Solar Azimuth Angle H5I_DATASET FLOAT
## 17 / Solar Zenith Angle H5I_DATASET FLOAT
## 18 / Surface Elevation H5I_DATASET FLOAT
## 19 / Visibility Index Map H5I_DATASET INTEGER
## 20 / Water Vapor Column H5I_DATASET INTEGER
## 21 / coordinate system string H5I_DATASET STRING
## 22 / flightAltitude H5I_DATASET FLOAT
## 23 / flightHeading H5I_DATASET FLOAT
## 24 / flightTime H5I_DATASET FLOAT
## 25 / fwhm H5I_DATASET FLOAT
## 26 / map info H5I_DATASET STRING
## 27 / to-sensor azimuth angle H5I_DATASET FLOAT
## 28 / to-sensor zenith angle H5I_DATASET FLOAT
## 29 / wavelength H5I_DATASET FLOAT
## dim
## 0 1
## 1 1
## 2 544 x 578 x 1
## 3 544 x 578 x 1
## 4 544 x 578 x 1
## 5 544 x 578 x 1
## 6 544 x 578 x 1
## 7 544 x 578 x 1
## 8 544 x 578 x 1
## 9 1
## 10 544 x 578 x 426
## 11 1
## 12 544 x 578 x 1
## 13 1
## 14 544 x 578 x 1
## 15 1
## 16 1 x 1
## 17 1 x 1
## 18 544 x 578 x 1
## 19 544 x 578 x 1
## 20 544 x 578 x 1
## 21 1
## 22 5732053
## 23 5732053
## 24 5732053
## 25 426 x 1
## 26 1
## 27 544 x 578 x 1
## 28 544 x 578 x 1
## 29 426 x 1
Import Data Dimensions
# open file for viewing-connection to file made (f for file)
fid <- H5Fopen(f)
fid
## HDF5 FILE
## name /
## filename
##
## name otype dclass
## 0 ATCOR_Input_File H5I_DATASET STRING
## 1 ATCOR_Processing_Log H5I_DATASET STRING
## 2 Aerosol Optical Depth H5I_DATASET INTEGER
## 3 Aspect H5I_DATASET FLOAT
## 4 Cast Shadow H5I_DATASET INTEGER
## 5 Dark Dense Vegetation Classification H5I_DATASET INTEGER
## 6 Haze-Cloud-Water Map H5I_DATASET INTEGER
## 7 Illumination Factor H5I_DATASET INTEGER
## 8 Path Length H5I_DATASET FLOAT
## 9 Processing Version H5I_DATASET STRING
## 10 Reflectance H5I_DATASET INTEGER
## 11 Shadow_Processing_Log H5I_DATASET STRING
## 12 Sky View Factor H5I_DATASET INTEGER
## 13 Skyview_Processing_Log H5I_DATASET STRING
## 14 Slope H5I_DATASET FLOAT
## 15 Slope_Aspect_Processing_Log H5I_DATASET STRING
## 16 Solar Azimuth Angle H5I_DATASET FLOAT
## 17 Solar Zenith Angle H5I_DATASET FLOAT
## 18 Surface Elevation H5I_DATASET FLOAT
## 19 Visibility Index Map H5I_DATASET INTEGER
## 20 Water Vapor Column H5I_DATASET INTEGER
## 21 coordinate system string H5I_DATASET STRING
## 22 flightAltitude H5I_DATASET FLOAT
## 23 flightHeading H5I_DATASET FLOAT
## 24 flightTime H5I_DATASET FLOAT
## 25 fwhm H5I_DATASET FLOAT
## 26 map info H5I_DATASET STRING
## 27 to-sensor azimuth angle H5I_DATASET FLOAT
## 28 to-sensor zenith angle H5I_DATASET FLOAT
## 29 wavelength H5I_DATASET FLOAT
## dim
## 0 1
## 1 1
## 2 544 x 578 x 1
## 3 544 x 578 x 1
## 4 544 x 578 x 1
## 5 544 x 578 x 1
## 6 544 x 578 x 1
## 7 544 x 578 x 1
## 8 544 x 578 x 1
## 9 1
## 10 544 x 578 x 426
## 11 1
## 12 544 x 578 x 1
## 13 1
## 14 544 x 578 x 1
## 15 1
## 16 1 x 1
## 17 1 x 1
## 18 544 x 578 x 1
## 19 544 x 578 x 1
## 20 544 x 578 x 1
## 21 1
## 22 5732053
## 23 5732053
## 24 5732053
## 25 426 x 1
## 26 1
## 27 544 x 578 x 1
## 28 544 x 578 x 1
## 29 426 x 1
# open the reflectance dataset-make connection to data (d for data)
did <- H5Dopen(fid,"Reflectance")
did
## HDF5 DATASET
## name /Reflectance
## filename
## type H5T_STD_I16LE
## rank 3
## size 544 x 578 x 426
## maxsize 544 x 578 x 426
# data is read as column,row,bands which is different from HDF5 Viewer
# grab dataset dimensions (x,y,z data) s for structure
sid <- H5Dget_space(did)
sid
## HDF5 DATASPACE
## rank 3
## size 544 x 578 x 426
## maxsize 544 x 578 x 426
# Import the dimesions of the data as column, row, band
dims <- H5Sget_simple_extent_dims(sid)$size
dims
## [1] 544 578 426
# close all open connections-otherwise could overwrite data
H5Sclose(sid)
H5Dclose(did)
H5Fclose(fid)
Read in Reflectance Data
# Extract slice of H5 file
# Import BAND 56 index is list of things to import (column,rows,bands)
b56 <- h5read(file = f,
name = "Reflectance",
index = list(1:dims[1],1:dims[2],56))
Convert Data to Matrix
# convert to matrix
b56 <- b56[,,1]
# plot data
image(b56)
image(log(b56),
main = "log transformed data")
# Look at distribution of data
hist(b56)
Clean the data
# assign no data value to object-will not be used in computations
b56[b56 == noDataValue] <- NA
# apply scale factor
b56 <- b56/scaleFactor
hist(b56)
Transpose Data (Rows and Columns)
# transpose the row and colums
b56 <- t(b56)
image(log(b56))
Create Spatial Extent
# Find data structure
str(mapInfo)
## chr [1(1d)] "UTM,1,1,325963.0,4103482.0,1.0000000000e+000,1.0000000000e+000,11,North,WGS-84,units=Meters"
## - attr(*, "Description")= chr "Basic Map information for envi style programs"
class(mapInfo)
## [1] "array"
# split out Map Info object
mapInfo <- strsplit(mapInfo, ",")
# Take out heirarchy strucutre in the data
mapInfo <- unlist(mapInfo)
# look at data structure
str(mapInfo)
## chr [1:11] "UTM" "1" "1" "325963.0" "4103482.0" ...
xMin <- as.numeric(mapInfo[4])
yMin <- as.numeric(mapInfo[5])
# to find the r coordinate for x and y max take min+(dimensions*resolution)
# to get spatial resolution
xres <- as.numeric(mapInfo[6])
yres <- as.numeric(mapInfo[7])
# bring in the coordinates to convert r coordinates to spatial coordinates
xMax <- xMin + (dims[1] * xres)
yMax <- yMin + (dims[2] * yres)
Create Extent Object
# Create extent (order xmin, xmax, ymin, ymax)
rasExt <- extent(xMin, xMax, yMin, yMax)
rasExt
## class : Extent
## xmin : 325963
## xmax : 326507
## ymin : 4103482
## ymax : 4104060
Create actual raster object
# raster with coordinate system extent
b56r <- raster(b56,
crs=CRS("+init=epsg:32611"))
# assign extent to raster
extent(b56r) <- rasExt
b56r
## class : RasterLayer
## dimensions : 578, 544, 314432 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 325963, 326507, 4103482, 4104060 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 0.0033, 0.5678 (min, max)
plot(b56r)
Import NEON functions
# install devtools
# install.packages("devtools")
library(devtools)
install_github("lwasser/neon-aop-package/neonAOP")
## Skipping install for github remote, the SHA1 (5b9a37fe) has not changed since last install.
## Use `force = TRUE` to force installation
library(neonAOP)
# use open band function to view data
b55 <- open_band(f,
bandNum = 55,
epsg = 32611)
b55
## class : RasterLayer
## dimensions : 578, 544, 314432 (nrow, ncol, ncell)
## resolution : 1, 1 (x, y)
## extent : 325963, 326507, 4102904, 4103482 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : in memory
## names : layer
## values : 0.0029, 0.5622 (min, max)
# plot data
plot(b55)
# import several bands (R band,G band,B band)
bands <-c(58, 34, 19)
epsg <- 32611
# create raster stack
RGBstack <- create_stack(f,
bands = bands,
epsg = epsg)
RGBstack
## class : RasterStack
## dimensions : 578, 544, 314432, 3 (nrow, ncol, ncell, nlayers)
## resolution : 1, 1 (x, y)
## extent : 325963, 326507, 4102904, 4103482 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:32611 +proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : X58, X34, X19
## min values : 0.0030, 0.0027, 0.0009
## max values : 0.5680, 0.4761, 0.3808
plot(RGBstack)
plotRGB(RGBstack,
stretch = "lin")
# cir image
bands <- c(90, 34, 19)
CIRstack <- create_stack(f,
bands = bands,
epsg = epsg)
plotRGB(CIRstack,
stretch = "lin")
